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Abstract. We analyze the synchronization transition of a globally coupled net¬ 
work of N phase oscillators with inertia (rotators) whose natural frequencies are 
unimodally or bimodally distributed. In the unimodal case, the system exhibits a 
discontinuous hysteretic transition from an incoherent to a partially synchronized 
(PS) state. For sufficiently large inertia, the system reveals the coexistence of a PS 
state and of a standing wave (SW) solution. In the bimodal case, the hysteretic syn¬ 
chronization transition involves several states. Namely, the system becomes coherent 
passing through traveling waves (TWs), SWs and finally arriving to a PS regime. 
The transition to the PS state from the SW occurs always at the same coupling, 
independently of the system size, while its value increases linearly with the inertia. 
On the other hand the critical coupling required to observe TWs and SWs increases 
with N suggesting that in the thermodynamic limit the transition from incoherence 
to PS will occur without any intermediate states. Finally a linear stability analysis 
reveals that the system is hysteretic not only at the level of macroscopic indicators, 
but also microscopically as verified by measuring the maximal Lyapunov exponent. 


1.1 Introduction 

The renowned Kuramoto model [14] for phase oscillators was generalized in 
1997 by Tanaka, Lichtenberg and Oishi (TLO) [24, 25] by including an ad¬ 
ditional inertial term. The TLO model revealed, at variance with the usual 
Kuramoto model, first order synchronization transitions even for unimodal 
distributions of the natural frequencies. TLO have been inspired in their ex¬ 
tension by a work of Ermentrout published in 1991 [7]; in this paper Er- 
mentrout has introduced a pulse coupled phase oscillator model with inertia 
to mimic the perfect synchrony achieved by a specific type of fireflies, the 
Pteroptix Malaccae (but also by certain species of crickets and humans). The 
peculiarity of these fireflies is that they are able to synchronize their flash¬ 
ing activity to some forcing frequency (even quite distinct from their own 
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intrinsic flashing frequency) with an almost zero phase lag. This happens be¬ 
cause they adapt their period of oscillation to that of the driving oscillator. 
After his introduction, the Kuramoto model with inertia has been employed 
to describe synchronization phenomena in crowd synchrony on Londons Mil¬ 
lennium bridge [23], as well as in Huygens pendulum clocks [4]. Furthermore, 
phase oscillators with inertia (rotators) have recently found application in the 
study the self-synchronization in power and smart grids [6,8,16,20,21], as 
well as in the analysis of disordered arrays of underdamped Josephson junc¬ 
tions [26]. Cluster explosive synchronization has been reported for an adaptive 
network of Kuramoto oscillators with inertia, where the natural frequency of 
each oscillator is assumed to be proportional to the degree of the corresponding 
node [12]. Rotators arranged in two symmetrically coupled populations have 
recently revealed the emergence of intermittent chaotic chimeras [17], imper¬ 
fect chimera states have been found in a ring with nonlocal coupling [11], and 
transient waves have been observed in regular lattices [13]. 

There is a wide literature devoted to coupled rotators with an unimodal 
frequency distribution, however only a really limited number of studies have 
been devoted to this model with a bimodal distribution, despite the subject 
being extremely relevant for the modelization of the power grids [8,18,20]. 
To our knowledge the synchronization transition in populations of globally 
coupled rotators with bimodal distribution has been previously analyzed only 
by Acebron et al. in [1]. More specifically, the authors considered a model with 
white noise and a distribution composed by two (5-functions localized at ±l?o- 
As suggested in [15], the presence of noise blurs the ^-functions in bell-shaped 
functions analogous to Gaussian distributions. Therefore one expects a similar 
phenomenology to the one observable for deterministic systems with bimodal 
Gaussian distributions of the frequencies. A multiscale analysis of the model, 
in the limit of sufficiently large j reveals the emergence from the incoherent 
state of stable standing wave solutions (SWs) and of unstable traveling wave 
solutions (TWs) via supercritical bifurcations, while partially synchronized 
stationary states (PSs) bifurcates subcritically from incoherence. However, the 
authors affirm that in the considered limit the bifurcation diagram coincides 
with that of the usual Kuramoto model without inertia [2]. 

In this article we analyze the synchronization transitions observable for 
unimodal and bimodal frequency distributions for a population of globally 
coupled rotators in a fully deterministic system. In particular, we will analyze 
the influence of inertia and of finite size effects on the synchronization transi¬ 
tions. Moreover, we will study the macroscopic and microscopic characteris¬ 
tics of the different regimes emerging during adiabatic increase and decrease 
of the coupling among the rotators. In particular, Section 1.2 will be devoted 
to the introduction of the model, of the indicators used to characterize the 
synchronization transition, and of the different protocols employed to perform 
adiabatic simulations. The Lyapunov linear stability analysis is introduced in 
Subsection 1.2.1. The results for unimodal distributions are reported in Sect. 
1.3, with a particular emphasis to the TLO mean field theory and its extension 
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to any generic state observable within the hysteretic region (SubSect. 1.3.2). 
The emergence of clusters of locked and whirling oscillators is described in 
details in SubSect. 1.3.3. The dynamics of the network for bimodal distribu¬ 
tions is analyzed in Section 1.4. In particular, SubSect. 1.4.1 is devoted to 
two non overlapping distributions and SubSect. 1.4.2 to largely overlapping 
Gaussian distributions. Sect. 1.5 report the result of linear stability analysis 
for the considered distributions. Finally, in Sect. 1.6 the reported results are 
briefly summarized and discussed. 


1.2 Model and Indicators 


By following Refs. [24,25], we study the following version of the Kuramoto 
model with inertia for N fully coupled rotators : 

m6^ + ei = + sm{ej - 6^) ( 1 . 1 ) 

3 


where 9i and f^i are, respectively, the instantaneous phase and the natu¬ 
ral frequency of the f-th oscillator, K is the coupling. In the following we 
will consider random natural frequencies 12^ Gaussian distributed accord¬ 


ing to: an unimodal distribution g{fi) — ^ with zero average and 

an unitary standard deviation or a bimodal symmetric distribution g{i2) = 

(f2-l7Q)2 (l7 + f2Q)2 

_ I ^-- 

2 ^/^ 


which is the overlap of two Gaussians with 


unitary standard deviation and with the peaks located at a distance 2l7o- 
To measure the level of coherence between the oscillators, we employ the 
complex order parameter [27] 






iOi 


( 1 . 2 ) 


where r{t) G [0 : 1] is the modulus and (j){t) the phase of the macroscopic 
indicator. An asynchronous state, in a finite network, is characterized by r ~ 
while for r = 1 the oscillators are fully synchronized and intermediate 
r-values correspond to partial synchronization. 

Another relevant indicator for the state of the rotator population is the 
number of locked oscillators Nl, characterized by a vanishingly small aver¬ 
age phase velocity Wj = and the maximal locking frequency 17 m, which 
corresponds to the maximal natural frequency jl7ij of the locked oscillators. 

In general we will perform sequences of simulations by varying adiabati- 
cally the coupling parameter K with two different protocols. Namely, for the 
first protocol (I) the series of simulations is initialized for the decoupled sys¬ 
tem by considering random initial conditions for {6i} and {w^}. Afterwards 
the coupling is increased in steps AK until a maximal coupling Km is reached. 
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For each value of K, apart the very first one, the simulations is initialized by 
employing the last configuration of the previous simulation in the sequence. 
For the second protocol (II), starting from the final coupling Km achieved by 
employing the protocol (I), the coupling is reduced in steps AK until K = 0 
is recovered. At each step the system is simulated for a transient time T/j fol¬ 
lowed by a period Tw during which the average value of the order parameter 
f and of the velocities {oJi}, as well as 17 m, are estimated. 

1.2.1 Lyapunov Analysis 

The stability of Eq. (1.1) can be analyzed by following the evolution of 
infinitesimal perturbations T = {SOi,... ,d9]\[,60i,... ,60^) in the tangent 
space, whose dynamics is ruled by the linearization of Eq. (1.1) as follows: 

K ■ . 

m S6i + ^ ^ cos (6j — 6i) {69j — 69i) . (1.3) 

^ i=i 

We will limit to estimate the maximal Lyapunov exponent Am, by employ¬ 
ing the method developed by Benettin et al. [3]. This amounts to follow the 
dynamical evolution of the orbit and of the tangent vector 'T for a time lapse 
Tw by performing Gram-Schmidt ortho-normalization at fixed time intervals 
At, after discarding an initial transient evolution Tfj. 

Eurthermore, the values of the components of the maximal Lyapunov vec¬ 
tor 'T can give important information about the oscillators that are more 
actively contributing to the chaotic dynamics. It is useful to introduce the 
following squared amplitude component of the normalized vector for each ro¬ 
tator [9,17] 

m , i = i,..., a . (i.4) 

The time average of this quantity gives a measure of the contribution of 
each oscillator to the chaotic dynamics. 


1.3 Unimodal Frequency Distribution 

1.3.1 Hysteretic Synchronization Transitions 

In Fig. 1.1 the results for a sequence of simulations obtained by following 
protocol (I) and (II) are reported for a not too small inertia (namely, m = 2) 
and unimodal frequency distribution. For the first protocol the system remains 
incoherent up to a critical value K = Kf ^ 2, where r jumps to a finite value 
and then increases with K reaching f ~ 1 for sufficiently large coupling. 
Starting from the last state and by reducing K one notices that f assumes 
larger values than during protocol (I) and the system becomes incoherent at a 
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Fig. 1.1. Unimodal frequency distribution, (a) Maximal locking frequency (blue 
triangles) and (b) time averaged order parameter f (black circles) as a function of 
the coupling K for two series of simulations performed following the protocol (I) 
(filled symbols) and the protocol (II) (empty symbols). The data refer to inertia: 
m = 2, for which we set AK = 0.2 and Km = 10; moreover N = 500, Tr = 5, 000 
and Tw = 200. The (magenta) diamonds indicate flp = ^for protocol (I) and 
the (green) squares fio = Kf for protocol (II). 


smaller coupling, namely itTf < Kf. This is a clear indication of the hysteretic 
nature of the synchronization transition in this case. 

For the chosen values of the inertia, we observe the creation of an unique 
cluster oi Nl locked oscillators with ~ 0, for larger m the things will 
be more complex, as we will discuss in the following. The maximal locking 
frequency f^M becomes finite for K > Kf and increases with K. The frequency 
attains a maximal value when f ~ 1, no more oscillators can be recruited 
in the large locked cluster. Once reached this value, even if K is reduced 
following the protocol (II), I7 m remains constant for a wide K interval. Then 
shows a rapid decrease towards zero by approaching Kf. This behavior 
will be explained in the following two sub-sections. 

1.3.2 Mean field theory 

In order to derive a mean field description of the dynamics of each single 
rotator, we can rewrite Eq. (I.l) by employing the order parameter definition 
(1.2) as follows 

mOi + 9i = fli — KrsiTL{6i — (p) ; (1-5) 

this corresponds to the evolution equation for a damped driven pendulum. 
Eq. (1.5) admits, for sufficiently small forcing frequency two fixed points: 

a stable node and a saddle. At larger frequencies f2i > flp the 
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saddle, via a homoclinic bifurcation, gives rise to a limit cycle. The stable 
limit cycle and a stable fixed point coexist until a saddle node bifurcation, 
taking place at fli = = Kr, leads to the disappearance of the fixed 

points and for 17^ > fijj only the limit cycle persists. This scenario is correct 
for sufficiently large inertia; at small m one has a direct transition from a 
stable node to a periodic oscillating orbit at Qi = fio = Kr [22]. Therefore 
for sufficiently large m there is a coexistence regime where, depending on the 
initial conditions, the single oscillator can rotate or stay quiet. The fixed point 
(limit cycle) solution corresponds to locked (drifting) rotators. 

The TLO theory [24, 25] has explained the origin of the first order hys- 
teretic transitions by considering two opposite initial states for the network: 
(I) the completely incoherent phase (r = 0) and (II) the completely synchro¬ 
nized one (r = 1). In case (I) the oscillators are all initially drifting with finite 
velocities Wi; by increasing K the oscillators with smaller natural frequencies 
\fii\ < Qp begin to lock {dJi = 0), while the other continue to drift. This is 
confirmed by the data reported in Fig. 1.1, where it is clear that the locking 
frequency f2M is well approximated by f2p. The process continues until all 
the oscillators are finally locked, leading to r = 1 and to a plateau in {2m- 

In the second case, initially all the oscillators are already locked, with an 
associated order parameter r = 1. Therefore, the oscillators can start to drift 
only when the stable fixed point solution will disappear, leaving the system 
only with the limit cycle solution. This happens, by decreasing K, whenever 
\f2i\ > f2p) = Kr. This is numerically verified, indeed, as shown in Fig. 1.1, 
where it is clear that the maximal locked frequency Hm remains constant 
until, by decreasing K, it encounters the curve flp and then f2M follows this 
latter curve down towards the asynchronous state. The case (II) corresponds 
to the situation observable for the usual Kuramoto model, where there is no 
bistability [14]. 

In both the considered cases there is a group of desynchronized oscillators 
and one of locked oscillators separated by a frequency, f2p {Qp ) in case (I) 
(case (II)). At variance with the usual Kuramoto model, both these groups 
contribute to the total level of synchronization, namely 


r = rL+rD 


( 1 . 6 ) 


where rp (rp) is the contribution of the locked (drifting) population. 

The contribution of the locked population is simply given by 

/ ^P,D 

cos^ 9g{Kr sin 9)d9 ; (1-7) 

-^P,D 

where 9p = sm~^{f2p/Kr) and 9p = sm~^{{2p/Kr) = 7r/2. 

The contribution rp oi the drifting rotators is negative and it has been 
estimated analytically by TLO by performing a perturbative expansion to 
the fourth order in l/(rnK) and l/(rnf2). The obtained expression, valid for 
sufficiently large inertia, reads as 
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(L8) 


with g{0) = g{—n). 

By considering an initially desynchronized (fully synchronized) system and 
by increasing (decreasing) K one can get a theoretical approximation for the 
level of synchronization in the system by employing the mean-field expression 

(1.7) , (1.8) and (1.6) for case I (II). In this way, two curves are obtained in 
the phase plane {K,r), namely (K) and r^^{K). For a certain coupling K 
the system can attain all the possible levels of synchronization between (K) 
and r^^{K). 

Let us notice that the expression for tl and ro reported in Eqs. (1.7) and 

(1.8) are the same for case (I) and (II), only the integration extrema change in 
the two cases. These are defined by the frequency which discriminates locked 
from drifting neuron, that in case (I) is I7p and in case (II) I7pi. It should be 
noticed that the value of these frequencies is a function of the order parameter 
r and of the coupling constant K, therefore one should solve implicit integrals 
to obtain r. 

However, one could also fix the discriminating frequency to some arbitrary 
value l7o and solve self-consistently the equations Eqs. (1.6), (1.7), and (1.8) 
for different values of the coupling K. This corresponds to solve the equation 

/ 9o roc 1 1 

cos^ 9g{Kr° sin 6)d9 — m / -— -^g{fl)dQ = — ; (1.9) 

-00 Joo K 


with 6*0 = sin~^{Q q/K r'^). A solution r° = r°(Ar, 12o) exists provided that 
l7o < ^d{K) = r^K. Therefore the part of the plane delimited by the curve 
will be filled with the curves r^i^K) obtained for different Qq values 
(as shown in Fig. 1.2(a)). These solutions represent clusters of iVp oscillators 
for which the maximal locking frequency and Nj^ do not vary upon changing 
the coupling strength. In particular, for K > ATf these states can be observed 
in numerical simulations in the portion of the phase space delimited by the 
two curves (K) and r^^{K) (see Fig. 1.2(b)). 


1.3.3 Clusters of Locked and Whirling Oscillators 

By observing the results reported in Fig. 1.2(a) for m = 6 , it is evident that 
the numerical data obtained by following the procedure (II) are quite well 
reproduced from the mean field approximation (solid grey curve). This is 
not the case for the theoretical estimation (dashed grey curve), which does 
not reproduce the step-wise structure revealed for the data corresponding to 
protocol (I). This step-wise structure emerges only for sufficiently large inertia 
(as it is clear from Fig. 1.1 (b), where it is absent for m = 2); this is due 
to the break down of the independence of the whirling oscillators: namely, 
to the formation of clusters of drifting oscillators moving coherently at the 
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Fig. 1.2. Unimodal Distribution. Panel (a): Average order parameter r versus the 
coupling constant K. Theoretical mean field estimates: the dashed (solid) grey curves 
refer to = »*£ +’'d = ’’’l! + ) 8-s obtained by employing Eqs. (1.7) and (1.8) 

following protocol (I) (protocol (11)); the (black) dot-dashed curves are the solutions 
r^{K, 17o) of Eq. (1.9) for different Qq values. The employed values from bottom to 
top are: Qq = 0.79, 1.09, 1.31 and 1.79. Numerical simulations: (black) filled circles 
have been obtained by following protocol (I) and then (II) starting from K = Q 
until Km = 20 with steps AK = 0.5; (black) empty triangles refer to simulations 
performed by starting from a final configuration obtained during protocol (I) and 
by decreasing the coupling from such initial configurations. The Panel (b) displays 
Nl vs K for the numerical simulations reported in (a). The numerical data refer to 
m = 6, N = 500, Tr = 5000, and Tw ~ 200. 


same non zero velocity [24]. Oscillators join in small groups to the locked 
stationary cluster and not individually as it happens for smaller inertia; this 
is clearly revealed by the behavior of versus the coupling K as reported 
in Fig. 1.2(b). 

Furthermore, once formed, these stationary locked clusters are particularly 
robust, as it can be appreciated by considering as initial condition a partially 
synchronized state obtained following protocol (I) for a certain coupling Ks > 
Ki. This state is characterized by a cluster of Nl locked; if now we reduce 
the coupling K, the number of locked oscillators remain constant until we 
do not reach the descending curve obtained with protocol (II), see the black 
empty triangles in Fig. 1.2. On the other hand r decreases slightly with K, this 
behavior is well reproduced by the mean field solutions of Eq. (1.9), namely 

r^{K, l7o) with l7o = ^p{Ks,r^(Ks)) = ^\j\ these are shown in Fig. 1.2 
(a) as black dot-dashed lines. As soon as, by decreasing K, the frequency 
l7o becomes equal or smaller than Qd, the order parameter has a rapid drop 
towards zero following the upper limit curve . These results indicate that 
hysteretic loops of any size are possible within the region delimited by the two 
curves f obtained by following protocol I and II respectively, as shown in [18]. 

For sufficiently small m, the synchronization occurs starting from the inco¬ 
herent state via the formation of an unique cluster of locked oscillators, and the 
size of this cluster increases with AT, as evident from Fig. 1.3 (a) for m = 2. 
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Fig. 1.3. Unimodal Distribution. Average phase velocity u>i of the rotators versus 
their natural frequencies Qi for N = 500 and inertia m = 2 (a) and m = 6 (b). In 
panel (a) (panel (b)) magenta triangles refer to K = 1.5 [K = 2.5), green diamonds 
to K = 2.5 {K = 5), red squares to K = 5.5 {K = 10) and black circles to 
K = 9.5 {K = 15). The insets report the time evolution of the order parameters 
r{t) for the corresponding coupling constants, apart for the extra blue line shown 
in the inset in (b) which refers to K = 1. For each simulation an initial transient 
Tr ~ 5, 500 has been discarded and the time averages have been estimated over a 
window Tw = 5, 000. 


At the same time the value of r also increases with K and its evolution is 
characterized only by finite size fluctuations vanishing in the thermodynamic 
limit (see the inset of Fig. 1.3 (a)). As already mentioned, the situation is 
different for sufficiently large inertia, now the partially synchronized phase is 
characterized by the coexistence of the main cluster of locked oscillators with 
dJi ~ 0, but also by the emergence of clusters composed by drifting oscillators 
with common finite velocities, see the data for dJi reported in Fig. 1.3 (b) 
for TO = 6. In particular, the clusters of whirling oscillators emerge always in 
couple and they are characterized by the same average velocity but opposite 
sign. These states are indicated as standing waves (SWs), therefore we have a 
SW coexisting with a partially synchronized stationary state (PS) (as shown 
in Fig. 1.3 (b). 

The effect of these extra clusters on the collective dynamics is to induce 
oscillations in the temporal evolution of the order parameter, as one can see 
from the inset of Fig. 1.3 (b). In presence of drifting clusters characterized 
by the same average velocity (in absolute value), as for to = 6 and K = 5 
in Fig. 1.3 (b), r exhibits almost regular oscillations and the period of these 
oscillations corresponds to the one associated to the oscillators in the drifting 
cluster. 
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1.4 Bimodal Distribution 

In this Section we consider a bimodal distribution, we will initially focus on 
two almost non overlapping Gaussians, namely we consider l7o = 2, while Sub- 
Section 1.4.2 is devoted to overlapping Gaussians, examined for I2o = 0.2. 

1.4.1 Non Overlapping Gaussians 

For l7o = 2 and sufficiently small inertia (m = 1 and 2), we observe a very rich 
synchronization transition, as shown in Fig. 1.4 (a). In particular, by following 
protocol (I) we observe that the system leaves the incoherent state abruptly by 
exhibiting a jump to a finite f value at ; above such value in the network 
emerges a single cluster of oscillators, drifting together with a finite velocity cs 
this corresponds to Traveling Wave (TW) solution . By further increasing 
K a second finite jump of the order parameter at denotes the passage 

to a Standing Wave (SW) solution, corresponding to two clusters of drifting 
oscillators with symmetric opposite velocities ±12o- A final jump at 
leads the system to a Partially Synchronized (PS) phase, characterized by an 
unique cluster of locked rotators with zero average velocity. By increasing the 
coupling the PS state smoothly approaches the fully synchronized regime . 
Starting from this final state the return sequence of simulations, following 
protocol (II), displays a simpler phenomenology. The network stays in the 
PS regime, characterized by an order parameter larger than that measured 
during protocol (I) simulations, until . For smaller coupling, 

the system leaves the PS state; however, depending on the realization of the 
natural frequencies and on the initial conditions, it can ends up or in a TW 
(most of the cases) or in a SW, or it can even reach directly the incoherent 
state (as shown in Fig. 1.4 (a) and Fig. 1.6 (a)). This first analysis clearly 
shows hysteretic effects and coexistence of macroscopic states with different 
level of synchronization for a wide range of couplings. 

Let us now try to examine the observed transitions in terms of the maximal 
locking frequency, I2 m- This frequency is now defined in a different way with 
respect to the unimodal distribution, in the present case Qm represents the 
maximal absolute value of the natural frequencies of the oscillators belonging 
to the main clusters present in the system, therefore in this estimation are 
considered both stationary and drifting clusters. As shown in Fig. 1.4 (b), 
Qm increases with K for protocol (I) simulations. In particular [2m shows a 
finite jump in correspondence oi K = , and then its evolution is rea¬ 

sonably well approximated by the curve I7p -|- f2o, where f2p = 4y/ 2^, By 

approaching the maximal frequency displays a constant plateau which 
extends beyond , this indicates that the two symmetric drifting clusters 
merge at K = giving rise to an unique locked cluster with zero average 
velocity, however no other oscillators join this cluster up to a larger coupling. 
Whenever this happens , Qm starts again to increase, but this time it follows 
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Fig. 1.4. Bimodal frequency distribution. Panel (a): Average order parameter f 
versus K for two series of simulations performed following the protocol (I) (filled 
symbols) and (II) (empty symbols). The dotted vertical blue line refers to 
the dashed-dotted magenta line to . Panel (b): Maximal locking frequency 
Qm (blue triangles) versus K for simulations reported in (a) for protocol (I) (filled 
symbols) and (II) (empty symbols). The magenta diamonds indicate Op = ^\/^ 
for protocol (I) and the (green) squares SId = Kf for protocol (II). The dashed 
magenta line represents the curve fip -I- Qo and the dashed green curve Qd -f I^o- 
In both panels the dashed orange vertical line denotes the critical value . The 
data refer to m = 2, Qq = 2, N = 2000, Tp = 10, 000 and Tw = 200, for the 
sequence of simulations we employed AK = 0.4 until K = 14.8 and AK = 0.4 for 
14.8 <K < 39.8. 


the curve Qp = “y Finally, for r ~ 1 the maximal locking frequency 
attains a maximal value. Moreover, by reducing the coupling, following now 
protocol (II), I7 m remains stacked to such a value for a wide K interval. The 
fully synchronized cluster is difficult to break down due to the inertia effects. 
Finally, Qm reveals a rapid decrease towards zero whenever it encounters the 
curve f2p) = Kf, initially it follows this curve, however as soon as the sys¬ 
tem desynchronizes towards a TW the decrease of 17 m is better described 
by the curve fip + 17o. The observed behavior can be explained by the fact 
that for K < {K < for protocol (I) (protocol (II)), the network 

behaves as two independent sub-networks each characterized by an unimodal 
frequency distribution, one centered at l7o and the other one at — l7o. The 
extension of the analysis reported in Section 1.3.2 for an unimodal distribu¬ 
tion not centered around zero simply amounts to shift the limiting curves I7p 
and l7p) by l7o. However, for sufficiently large coupling constant, once the 
system exhibits only one large cluster with zero velocity, the network behaves 
as a single entity and 17 m closely follows I7p or flp as for a single unimodal 
distribution centered in zero. 

Let us now describe the TW and SW states in more details with the help 
of the examples reported in Fig. 1.5 for inertia m = 2 and N = 2000. The 
TW is an asymmetrical cluster of whirling oscillators with a finite velocity 
dJi ~ l7o, in particular in Fig. 1.5 (a) the oscillators have natural frequencies 
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Fig. 1.5. Bimodal frequency distributions. Average phase velocity txn of the rotators 
versus their natural frequencies fli for coupling strength K=6.2 (a) and K=6.8 
(c). Panels (b) and (d) display the order parameter r(t) (black line) versus time 
for the same coupling constants as in (a) and (c), respectively. The dashed black 
(continuous grey) line denotes the time evolution of {vn). For each simulation an 
initial transient time Ta = 1000 has been discarded and the average are estimated 
over a time interval Tw ~ 200. In both cases m = 2, N = 2000 and Oq — 2. 


in a range around i7o, namely 0.67 < f2i < 3.34. The effect of this cluster on 
the collective dynamics is to increase the average value of the order parameter 
without inducing any clear oscillating behavior in r{t). However, oscillators 
with positive natural frequencies are much more synchronized with respect to 
the ones with negative frequencies, as can be inferred by observing the order 
parameter Vp (jn) estimated only on the the sub-population of oscillators 
with positive (negative) natural frequencies and reported in Fig. 1.5 (b). We 
believe that the emergence of the TW state is related to the finite sampling 
of the distribution of the natural frequencies, which due to finite size effects 
can be non perfectly symmetric. The asymmetric cluster will emerge around 
+^0 {—^2o) depending on the positive (negative) sign of the average natural 
frequency. 

The SWs are observable at larger coupling constants, this is character¬ 
ized by two symmetrical clusters with opposite average velocities ~ ±I2o) as 
shown in Fig. 1.5 (c). The presence of these two clusters now induces clear 
periodic oscillations in the order parameter r(t), as observable in Fig. 1.5 (d). 
The period of the oscillations is related to |l7ol; i-e. the average frequency 
of the clustered oscillators. However, at variance with the results reported 
in Fig. 1.3 (b) for the unimodal distribution the two symmetric clusters do 
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not coexist with a cluster of locked oscillators with zero average velocity. By 
examining separately and reported in Fig. 1.3 (d), we notice that each 
sub-population is much more synchronized than the global one, in fact Vp and 
Vn have a higher average value than r with superimposed irregular oscillations. 



Fig. 1.6. Bimodal frequency distribution with f?o = 2. Average order parameter r 
versus K for various system sizes N-. (a) m = 1, (b) m = 6. The numerical data 
have been obtained by following protocol (I) and then protocol (II) from K — 0 
up to Km = 20 {Km = 200) for inertia m=l (m=6) with AK — 0.2 {AK = 0.5). 
The vertical dashed blue line refers to K^^. The average has been performed over 
a time window Tw = 200, after discarding a transient time Tr = 5, 000 — 50, 000 
depending on the system size; the larger Tr have been employed for the larger N. 


The data reported so far refer to a single system size, however finite size 
effects are quite relevant for this model, as shown in [18] for unimodal distri¬ 
butions. In Fig. 1.6 (a) we report the synchronization transition for several 
system sizes, namely 1,000 < N < 50,000 for a small inertia value (m = 1). 
We observe that and increase with the size iV; in particular, the 

incoherent state is observable on a wider coupling interval by increasing N 
(similarly to what reported in [18] for unimodal distributions). Finite size fluc¬ 
tuations induce transitions from the incoherent branch to the TW branch and 
from this to the SW branch. The fact that we do not observe transition back 
to the original states indicates that the energy barriers are higher from these 
sides. A quite astonishing result is the fact that the transition value and 
seem completely independent from N. The combination of these results 
seem to suggest that in the thermodynamic limit the incoherent state will 
loose stability at and therefore the two branches corresponding to TW 
and SW will be no more visited, at least by following protocol (I). 

By observing all the data reported in Fig. 1.6 (a) for various N and for 
protocol (I) and (II), it seems that there are clear indications that the two 
branches of solutions, corresponding to TW and SW, emerge via a supercrit¬ 
ical bifurcation at the same coupling, namely K ~ 3.8, while the transition 
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60 h 
PS 


Fig. 1.7. Critical value as a function of the inertia. The dashed line represents 
the fit of the numerical data and indicates a linear increasing of as a function of 
the inertia being the fit = 4(1.245 + 1.0525m). For all cases N = 2000, Qq = 2. 
The data have been obtained by employing protocol (I) and for each simulation an 
initial transient time Tr = 5000 has been discarded and data are averaged over a 
time Tw = 200. 


to PS is clearly subcritical. These results confirm the analysis reported in [1] 
for a system with noise (in particular, see Fig. 17 in that paper). However, 
Acebron et al. affirm that the SW is stable, while the TW is unstable. From 
our results, both branches seem to become inaccessible (in absence of noise) 
from the incoherent state, while at least a part of these branches appear to be 
reachable from the PS state by decreasing K below following protocol 

(II). Another important difference with respect to the results reported in [1] 
is that the PS regime is clearly hysteretic revealing two coexisting branches 
of PS states visited by following protocol (I) or (II). 

As shown in Fig. 1.6 (b), for larger inertia (namely, m = 6), the transition 
from the incoherent state following protocol (I) occurs via the emergence of 
many small clusters leading finally to a SW state. In this case the critical value 
at which the incoherent state looses stability seems to saturate to a constant 
value already for Af > 2, 000. The value of is also in this case insensible 
to the system size. For large inertia values, the TWs seem no more observable. 

As a further aspect, we will report the numerical results of the dependence 
on the inertia of the critical coupling constant , while the value of cs 
4.9 is independent not only by N, but also by the inertia. As shown in Fig. 
1.7, increases linearly with the inertia and this scaling is already valid 
for not too large inertia values. The linear scaling with the inertia is analogous 
to the scaling recently found within a theoretical mean-field analysis for the 
coupling , which delimits the range of linear stability of the asynchronous 
state [I, 10]. In particular, the authors in [18] have shown for a Gaussian 
unimodal distribution of width a that ~ 2 ct(0.64 + toct), which shows a 

linear dependence on the inertia and a quadratic dependence on the variance 
of the frequency distribution. 
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In the final part of this sub-section we perform an analysis analogous to 
that reported in Sub-Sect. 1.3.3, in particular starting from states with a 
finite level of synchronization obtained by following protocol (I) we decrease 
the coupling and observe how these states evolve. In Fig. 1.8, we report the 
results of these simulations (shown as empty triangles) for two different inertia 
values, namely m = 1 and m = 10. Starting from PS states we observe that 
the cluster survives until the descending curve obtained with protocol (II) is 
encountered, analogously to the results reported in Fig. 1.3 for the unimodal 
distribution. Therefore any part of the hysteretic portion of the {K, r)-plane 
delimited by the PS curves obtained via protocol (I) or (II) is accessible . 
However, if one starts for m = 1 from a TW or a SW state, one observes only 
two curves (corresponding to the TW and SW branches previously discussed) 
which seem to end up at the same critical coupling which is smaller than 
. Therefore it seems that there are no evidences of hysteresis for this 
small inertia for SW and TW solutions (as shown in Fig. 1.8 (a)). For large 
inertia values m = 10, since now, apart the SW solutions, there are solutions 
with many small clusters, the situations is much more complex. By starting 
from different values oi K < and by decreasing K, these curves seem all 
to end up at the same critical coupling smaller , see Fig. 1.8 (b). These 
results suggest that for large inertia values is possible to observe a continuum 
of possible states even starting from states characterized by (many) drifting 
clusters and that these states coexist in a wide range of coupling. 




Fig. 1.8. Bimodal frequency distribution. Average order parameter f versus the 
coupling constant K for m = 1 and N = 10,000 (panel (a)) and for m = 10 and 
N = 2, 000 (panel (b)). The filled circles have been obtained by following protocol (I) 
and then (II) starting from K = 0 until Km = 20 {Km = 100) with steps AK = 0.2 
(AK = 0.5); the empty triangles refer to simulations performed by starting from 
a final configuration obtained during protocol (I) and by decreasing the coupling 
from such initial configurations. The numerical data refer to f7o = 2, Ta = 50000 
{Tr = 5000), and Tw = 2000. 
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1.4.2 Overlapping Gaussians 

In this sub-section, we analyze a bimodal distribution, where the two Gaus¬ 
sians are largely overlapping, since l7o = 0.2. In this case we expect to observe 
a phenomenology of the synchronization transition quite similar to the one 
seen for the unimodal case. In Fig. 1.9(a) is reported the average order pa¬ 
rameter f versus the coupling constant K estimated by following protocol 
(I) and (II) for various inertia values and for a fixed system size, namely 
N = 10,000. We observe that all the curves obtained for protocol (II) almost 
overlap irrespectively of the used inertia, while the protocol (I) curves reveal 
a strong dependence on m. In particular, the hysteretic region widens with 
TO. For small inertia values, namely to = 1 and 2, there is a sudden transition 
from the asynchronous state to a PS state at and neither traveling waves 
nor standing waves are observable: a single cluster at zero velocity emerges 
in correspondence of and the order parameter never shows oscillating 
behavior in time. 

For TO = 6, it is possible to observe a scenario similar to the one reported 
in Fig. 1.3 (b), where not only a cluster at zero velocity is present, but also 
two symmetrical clusters at finite velocities emerge. In particular, following 
protocol (I) for K > 2.4 a small cluster of locked oscillators emerges; at 
larger coupling, namely K > 3, two symmetrical clusters of whirling oscillators 
emerge and coexist with the zero velocity cluster. Finally, at K = 10.8 the 
PS regime arises, corresponding to a single large cluster of locked oscillators. 
Furthermore, in the range Z < K < 10.8 the order parameter reveals irregular 
oscillations. A more detailed analysis is needed to understand the origin of 
these oscillations as done in the next Section. 

Finally, we examine the influence of the system size on the studied tran¬ 
sitions: for TO = 1 the results for the protocol (I) [protocol (II)] simulations 
are reported in Fig. 1.9(b) for sizes ranging from N = 1, 000 to TV = 10, 000. 
It is immediately evident that the transition from synchronized state to the 
asynchronous state, following protocol (II), does not depend on N: for all con¬ 
sidered sizes the transition happens in correspondence of AT ~ 2, analogously 
to what reported for unimodal distributions [18]. Starting from the incoher¬ 
ent regime and following protocol (I) the system reveals a jump to a finite r 
value for critical couplings increasing with fV, quite similar once more to the 
results reported for unimodal distributions. We can conclude this sub-section 
by afhrming that the phenomenology seen for bimodal, but largely overlap¬ 
ping, distributions should not differ much from the one observed for unimodal 
distributions. 


1.5 Linear Stability Analysis 

To better characterize the synchronization transitions and the stability of the 
observed states it is worth estimating the maximal Lyapunov exponent Am 
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Fig. 1.9. Bimodal frequency distribution for i?o = 0.2. Panel (a) : Average order 
parameter f versus K for various inertia values and N=1000. The numerical data 
have been obtained by following protocol (I) and then protocol (II) from K = 0 
up to Km = 20 for all inertia values with AK = 0.2 . Panel (b): Average order 
parameter r versus the coupling constant K for various system sizes N and m=l. 
Data have been obtained by averaging the order parameter over a time window 
Tw ~ 200, after discarding a transient time Tr = 5,000 — 50, 000 depending on the 
system size. 


following protocol (I) and (II) for an unimodal and a bimodal distributions. 
This quite time consuming analysis has been performed for a single inertia 
value m = 6 and a single system size = 1, 000, the scaling of Am with N 
will be discussed in the following for specific coupling constant values. 

In general, we observe that once the system fully synchronizes, Am van¬ 
ishes; therefore for most of the simulations associated to protocol (II) corre¬ 
sponding to fully synchronized cluster down to the desynchronization transi¬ 
tion, Am is zero. This is not the case for protocol (I) simulations which reveal 
a positive Am as soon as f is non zero. Thus indicating that not only the 
dynamics characterized in terms of the macroscopic order parameter f is hys- 
teretic, but also at the level of the microscopic dynamics, investigated via Am, 
the system has a clear hysteretic behavior. 

The behavior of Am with K exhibits chaotic dynamical states with win¬ 
dows of regularity for both unimodal and bimodal distribution with 17 = 2, as 
shown in Figs. 1.10 (a) and (b). As a general aspect, we observe the maximal 
level of chaoticity immediately after the transition from the incoherent state 
to partially coherence, where small clusters of synchronized oscillators and 
drifting oscillators coexist. The increase of r is accompanied by a trend of Am 
to decrease and finally to vanish for f —>■ 1 

An important aspect to understand is if this dynamics is weakly chaotic 
or not, in particular this amounts to verify if, in the thermodynamic limit. 
Am will vanish or will remain finite. In order to test for this aspect, we have 
considered a configuration obtained by following protocol (I) for a specific 
coupling and analyzed Am versus the system size for 200 < N < 32,000. The 
results for unimodal distributions, as well as for bimodal ones with l7o = 2 
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Fig. 1.10. Maximal Lyapunov exponent Am and the average order parameter r 
versus K for unimodal (bimodal with 17o = 2) are shown in panel (a) (panel (b)) 
and in panel (c) (panel (d)), respectively. The numerical data have been obtained 
by following protocol (I) (black circles) and then protocol (II) (red diamonds) from 
A” = 0 up to Km = 30 {Km = 80) with AK = 0.2 {AK — 0.8). For panels (a),(c) 
Tr = 500, and Tw = 50,000; for panels (b),(d) Tr = 500, and Tw = 400000. 
The different symbols in (a) and (b) denote the value for which the further analysis 
reported in panel (e) has been done. Panel (e): Am versus N for different couplings 
and frequency distributions. Blue circles refer to unimodal distribution and coupling 
constant K — 6.5; magenta squares (green diamonds) refer to binomial distributions 
with = 0.2 and K — 6.7 (17o = 2 and K = 9.5). Am has been averaged over a time 
window Tw = 4,000 — 400, 000, after discarding a transient time Tr = 1,000 — 10, 000 
depending on the system size. For all panels m = 6 and N = 1, 000. 


and = 0.2 are shown in Fig. 1.10 (e). It is clear for all the considered cases 
that the system remains chaotic for diverging system sizes. 

As a final aspect we would like to understand which oscillators contribute 
more to the chaotic dynamics of the system; this can be understood by mea¬ 
suring the average squared amplitude of the components of the maximal Lya¬ 
punov vector (see the definition reported in Eq. 1.4). In particular, we 
consider the three cases analyzed in Fig. 1.10 (e) for N = 1,000. The cor¬ 
responding results are shown in Fig. 1.11. From panel (a) and (b) of the 
figure it is clear that for the unimodal distribution, as well as for the largely 
overlapping bimodal distributions, the chaotic activity is associated almost 
exclusively to the rotators which are outside the large clusters of locked oscil- 
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a. 

Fig. 1.11. Average squared amplitude of the components of the maximal Lyapunov 
vector (black circles) and average frequencies oJi (grey diamonds) versus the ro¬ 
tator index for unimodal (a) and bimodal distribution with Oq — 0.2 (b). In both 
panels the oscillators are ordered according to the values of Panel (c); (black 
circles) and average frequencies oJi (grey diamonds) of the oscillators as a function 
of their natural frequency for bimodal distribution with f2o = 2. Panel (a) refers to 
coupling constant K = 6.5, panel (b) to if = 6.7 and panel (c) to if = 9.5. For all 
panels inertia m = 6 and N = 1, 000, Tr = 500, and Tw ~ 400000. 


lators with Wi ~ 0. Thus confirming recent results reported for two coupled 
populations of rotators with identical natural frequencies [17]. 

However, the situation for the bimodal distribution with l7o = 2 is dif¬ 
ferent; in particular, as shown in Fig. 1.11 (c), the network for this large 
value of the inertia and the considered coupling does not exhibit a cluster 
of locked oscillators with u)i ~ 0, but only drifting clusters. In this case the 
rotators outside and inside the clusters seem to contribute equally to the max¬ 
imal Lyapunov vector, with the possible exclusion of a group of rotators with 

— 12q. 

1.6 Conclusions 

We have studied the synchronization transition for a globally coupled Ku- 
ramoto model with inertia for different frequency distributions. For the uni¬ 
modal frequency distribution we have shown that clusters of locked oscillators 
of any size coexist within the hysteretic region. This region is delimited by 
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two curves in the plane individuated by the coupling and the average value of 
the order parameter. Each curve corresponds to the synchronization (desyn¬ 
chronization) profile obtained starting from the fully desynchronized (syn¬ 
chronized) state. For sufficiently large inertia values, clusters composed by 
drifting oscillators with opposite velocities (standing wave state) emerge in 
addition to the locked oscillators clusters. The presence of clusters of whirling 
rotators induces oscillatory behavior in the order parameter. 

For bimodal frequency distribution the scenario can become more complex 
since it is possible to play with an extra parameter: the distance between the 
peaks of the distributions. For simplicity we have analyzed only two cases: 
largely overlapping distributions (f2o = 0.2), and almost not overlapping dis¬ 
tributions (f2o = 2). The phenomenology observed for Hq = 0.2 resembles 
strongly that found for the unimodal distribution. The analysis of the non 
overlapping case reveals new interesting features. In particular, the transi¬ 
tion from incoherence to coherence occurs via several states: namely, travel¬ 
ing waves, standing waves and finally partial synchronization. This scenario 
resembles that reported for the usual Kuramoto model for a bimodal distribu¬ 
tion [5,15,19]. However, in our case the transition is always largely hysteretic, 
and for non overlapping distributions, traveling waves are clearly observable 
at variance, not only with the results for the Kuramoto model [5,15,19], but 
also with the theoretical phase diagram reported in [1] for oscillators with 
inertia. A peculiar aspect is that in the thermodynamic limit we expect a 
direct discontinuous jump from the incoherent to the coherent phase, without 
passing through any intermediate state. The critical coupling required 

to pass from incoherence to partial synchronization is independent of the sys¬ 
tem size and grows linearly with inertia, while the partially synchronized state 
looses its stability at a smaller coupling which is the same for 

any inertia value and system size. 

Finally, by performing a linear stability analysis we have been able to show 
that the hysteretic behavior is not limited to macroscopic observables, as the 
level of synchronization, but it is revealed also by microscopic indicators as the 
maximal Lyapunov exponent. In particular, we expect that in a large interval 
of coupling values chaotic and non chaotic states will coexist. 
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